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Abstract. We study the structure of neutron stars in perturbative f{R) gravity models with 
realistic equations of state. We obtain mass-radius relations in a gravity model of the form f{R) = 
R+aE?. We find that deviations from the results of general relativity, comparable to the variations 
due to using different equations of state (EoS'), are induced for \a\ ~ cm^. Some of the soft 
EoS' that are excluded within the framework of general relativity can be reconciled with the 2 solar 
mass neutron star recently observed for certain values of a within this range. For some of the EoS' 
we find that a new solution branch, which allows highly massive neutron stars, exists for values 
of a greater than a few 10^ cm^. We find constraints on a for a variety of EoS' using the recent 
observational constraints on the mass-radius relation. These are all 5 orders of magnitude smaller 
than the recent constraint obtained via Gravity Probe B for this gravity model. The associated 
length scale ^/a ~ 10^ cm is only an order of magnitude smaller than the typical radius of a neutron 
star, the probe used in this test. This implies that real deviations from general relativity can be 
even smaller. 
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1 Introduction 

The current accelerated expansion of the universe has been confirmed by many independent ob- 
servations. The supporting evidence comes from the supernovae la data [1-3], cosmic microwave 
background radiation [4-6], and the large scale structure of the universe [7, 8]. Although the cos- 
mological constant is arguably the simplest explanation and the best fit to all observational data, 
its theoretical value predicted by quantum field theory is many orders of magnitude greater than 
the value to explain the current acceleration of the universe. This problematic nature of cosmo- 
logical constant has motivated an intense research for alternative explanations and the reasonable 
approaches in this direction can be divided into two main categories, both of them introducing new 
degrees of freedom [9]: The first approach is to add some unknown energy- momentum component 
to the right hand side of Einstein's equations with an equation of state p/p ^ —1, dubbed dark 
energy. In the more radical second approach, the idea is to modify the left hand side of Einstein's 
equations, so-called modified gravity. Trying to explain such perplexing observations by modifying 
gravity rather than postulating an unknown dark energy has been an active research area in the 
last few years and in this paper we adopt this path. 

A modified theory of gravity has to explain the late time cosmology, and also be compatible 
with the constraints obtained from solar system and laboratory tests. However, it is not easy to 
construct theories of gravity with these requirements. Nevertheless, a class of theories, called f{R) 
models [10-12], has attracted serious attention possibly because of its (deceptive) simplicity. Today 
there exist viable f{R) models which are constructed carefully to be free of instabilities, and to 
pass the current solar system and laboratory tests [13-20]. 

The strong gravity regime [21] of these theories is another way of checking their viability. 
In this regime, divergences stemming from the functional form of f(R) may prevent the existence 
of relativistic stars in such theories [22-28], but thanks to the chameleon mechanism the possible 
problems jeopardizing the existence of these objects may be avoided [29, 30]. Furthermore, there 
are also numerical solutions corresponding to static star configurations with a strong gravitational 
field [31, 32] where the choice of the equation of state for the star is crucial for the existence 
of solutions, and therefore a polytropic equation of state is used in these works to overcome the 
possible problems related to the equation of state. 

Another approach to probe the viability of f(R) theories in the strong gravity regime is to use 
a method called perturbative constraints, or order reduction [33, 34]. The motivation behind using 
this technique in the present context is the thought that the reason of all the problems encountered 
in modifying gravity may be the outcome of considering these modifications as exact ones. The 



- 1 - 



main issue with the exact modifications are the problems arising in curvature scales which are not 
originally aimed by these modifications. In the perturbative constraints approach, the modifications 
are viewed as next to leading order terms to the terms coming from Einstein's General Relativity. 
Treating f(R) gravity via perturbative constraints at cosmological scales is considered in [35, 36] 
and for compact objects in [37]. In this manuscript we further examine the existence and properties 
of relativistic stars in the context of f(R) models via perturbative constraints. 

Specifically, we study a f{R) model of the form f{R) = R + aE? and constrain the value of a 
with the recent constraints on the mass-radius relation [38]. Such a gravity model in the weak field 
limit is known to reduce to Yukawa-like potentials and has been recently constrained by binary 
pulsar data [42] as a < 5 x 10^^ cm^. Here we find that in the strong gravity regime the constraint 
on perturbative parameter is a < 10^*^ cm^. This value does not contrast with the value obtained 
in [42] as they argue that the value of a could be different at different length scales. 

The plan of the paper is as follows: In section II, we assume a perturbative form of f{R) 
modified gravity model and obtain the field equations. Assuming also perturbative forms of metric 
functions and the hydrodynamical quantities we obtain the modified Tolman-Oppenheimer-VolkofF 
(TOV) equations. In section III, the modified TOV equations are solved numerically for various 
forms of the equation of state, the functional form of /(R), and various values of the perturbation 
parameter a. Finally, in the discussion section, we comment on the results of numerical study and 
on the significance of the scale of the perturbation parameter a. 

2 Modified TOV equations of f{R) gravity 

The action of f{R) gravity models is the simplest generalization of the Einstein-Hilbert action. 
Here, instead of having a linear function of Ricci scalar as the Lagrangian density we have a 
function of it: 



where g denotes the determinant of the metric g^j^, and R is the Ricci scalar. We set G = 1 and 
c = 1 in the action and in the rest of this section. Here we are considering the metric formalism of 
gravity and therefore matter only couples to the metric, and the Levi-Civita connection is derived 
from the metric. 

A straightforward variation of the action (2.1) with respect to the metric gives fourth order 
differential equations of g^^- However, treating higher than second order differential equations in 
4 dimensions is problematic. For this reason, we adopt a perturbative approach as suggested in 
[37] and choose f{R) such that all terms higher than second order will be multiplied by a small 
parameter a. The meaning of smallness of the parameter a is explained in the next section. 

Since this is a perturbative approach, the action and the field equations must have the form 
of those of general relativity for a = 0. So we choose the function f{R) to have the form 



without a constant piece, i.e. without a cosmological constant. Here h{R) is, for now, an arbitrary 
function of R and 0{a^) denotes the possible higher order corrections in a. Variation of the action 
(2.1) with respect to the metric, with the form of f{R) as given in (2.2), results in field equations 
which are 




(2.1) 




(2.2) 



(1 + ahR)G, 



-a{h - hRR)gi 



(2.3) 



- 2 - 



where G^y = R^i, — ^Rg^u is the Einstein tensor and = ^ is the derivative of h(R) with respect 
to the Ricci scalar. 

We are interested in sphericahy symmetric solutions of these field equations inside a neutron 
star, so we choose a spherically symmetric metric with two unknown, independent functions of r, 

ds^ = -e^'*'"dt^ + e^^'^dr'^ + r'^{de'^ + sin^ ed(j?). (2.4) 

Perturbative solution of the field equations means that g^y can be expanded perturbatively in 
a and therefore the metric functions have the expansions (pa = + «(/>! + • • • and Aq, = A + aAi + . . .. 
The energy-momentum tensor, which is on the right hand side of the field equations, is the energy- 
momentum tensor of the perfect fiuid. For a perturbative solution, then, hydrodynamics quantities 
are also defined perturbatively: Pa = P + oipi + . . . and Pa = P + olPi + . . .. Note that we denote 
zeroth order quantities, which can be obtained solving Einstein equations, without a subscript. 

Then the "tt" and "rr" components of field equations become 

- Si^pa = + e-2-^'^(l - 2rA'^)r-2 + ahRi-r'"^ + 6-^^(1 - 2rX')r-'^) 

-^a{h - HrR) + e-2^a[/i'^r-i(2 - rA') + h'^l (2.5) 

87rP„ = -r"^ + e-2-^'^(l + 2r0^)r-2 + ahR{-r-'^ + 6-^^(1 + 2r0')r"^) 

-^a(/i - HrR) + e-'^^ah'Rr-^{2 + rcj)'), (2.6) 

where prime denotes derivative with respect to radial distance, r. 

We would like to solve pa and Pa up to order a. On the right hand side, terms containing 
hR and its derivatives are already first order in a. Therefore, for the quantities that are multiplied 
by them, we use the zeroth order quantities which can be obtained from the Einstein equations 
written for this metric. In the equations above, hence, in terms that are multiplied by a we have 
already written zeroth order quantities. Using the "tt" and "rr" components of Einstein equations, 

-Svrp = + 6-2^(1 -2rA')r-^ (2.7) 

SvrP = -r~2 + e~2^(l + 2r0>-^ (2.8) 

we replace terms multiplied by Hr in (2.5) and (2.6) with —Sirp and SvrP, respectively. Additionally, 
combining a order terms on the left hand sides of (2.5) and (2.6) we obtain. 



2" - 1 - e-^-^'^d - 2rA' 



57rr Pa 

2 



+ahRr 



(2.9) 



Bvrr^P^ = -1 + e-2^"(l + 2r(l)'^) 

(2.10) 



+ahRr'^ 



SvrP - i f A _ + e-2V-i(2 + rcl>')^ 
n-R J hR 



To define a mass parameter, we assume a solution that has the same form of the exterior 
solution for the metric function Aq,. This form of solution has been previously suggested by the 
authors of [37]. Therefore we define 

')\ Ma , , 

e"2^" = 1 -. (2.11 

r 

Here, similar to pa, Ma is expanded in a as Ma = M + aMi + . . ., where M is the zeroth order 
solution, which, in general relativity, is given in terms of p as 

M = 8ir p{r)r^dr. (2.12) 
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Taking the derivative of Ma with respect to r one obtains 



dr 



1 



-2Ac 



{l-2r\'^). 



Substituting this into (2.9) and arranging terms, one gets the first modified TOV equation, 
dMa 



dr 



iirr'^Pa 



ahR 



+ (4^pr3-2r + fM)^ 



(2.13) 



(2.14) 



To obtain this equation, we also substitute general relativistic form of (2.11), e""^^ = 1 — 
The conservation equation of energy-momentum tensor of a perfect fluid, V^T^j^ 
equivalent to the hydrostatic equilibrium equation. 



M 
r 



0, is 



dPg 

dr 



-{pa + Pa 



dr 



(2.15) 



Therefore in order to get the second modified TOV equation we pull from the "rr" field equation, 
eq.(2.10). After some straightforward manipulations one gets. 



2(r-M„): 



dr 



2 Ma 

Svrr Pa H ahji 

r 



8^r'P+^4{^-R) 
+ (2r - |M + 47rPr3)J| 



(2.16) 



Note that when one sets a to zero, one gets the original TOV equations for general relativistic 
quantities. Similar to the case in general relativity, the modified TOV equations, (2.14), (2.15) and 
(2.16), are solved numerically for some special functional form of h{R). Obviously perturbation 
expansion parameter a introduces a new scale into the theory. By choosing a realistic equation 
of state we compute mass-radius relation for various values of a and this way we put a bound on 
a for perturbative f{R) gravity models with various forms of h{R).* This numerical analysis is 
explained in the next section. 



3 Numerical model and astrophysical constraints on the value of a 

Eqs. (2.14) and (2.16) describe any spherical mass distribution in a general f{R) theory in pertur- 
bative approach. Interesting results would be obtained in the case of neutron stars which have the 
highest compactness ratio r] = 2GM*/c^i?* and curvatures ^ = GM^./c^Rl [21]. In order to spe- 
cialize these equations for describing neutron stars they must be supplemented by an appropriate 
equation of state (EoS). However, the EoS of nuclear matter at the densities prevailing in neutron 
stars is not very well constrained by nuclear scattering data and there is a number of EoS leading 
to different mass-radius (M-R) relations for neutron stars. 

For integrating the TOV equations in GR, it is possible to interpolate the tabulated relation 
between density and pressure. In f{R) theories, where one needs high order derivatives of pressure 
with respect to density, interpolation leads to numerical problems. In order to circumvent this 
problem, we employ analytical expressions obtained by fitting the tabulated data. Such analytical 
representations are already provided by [53] for two EoS', EPS and SLy. For the rest of the EoS', 
namely AP4, GSl, MPAl, and MSI, we used analytical representations provided by [54] obtained 
by fitting the tabulated data with a function which is an extension of the function provided in 

*A realistic equation of state employs different physics for different regions of the neutron star, as opposed to a 
single polytropic relation prevailing throughout the star. 
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[53]. The 6 EoS' we have chosen constitute a representative sample (see figure 1 in [38]) for muhi- 
nucleonic and condensate inner composition possibilities (see [40] for the description of all these 
EoS'). We have not employed any strange quark matter EoS. As such stars are not gravitationally 
bound, alternative gravity models does not produce different mass-radius relations for such objects. 

We numerically integrate Eqs. (2.14), (2.15) and (2.16) supplemented by the analytical ex- 
pression for the EoS, employing a Runge-Kutta scheme with fixed step size of Ar = 0.001 km. We 
obtain a sequence of equilibrium configurations by varying the central density pc from 2 x lO'^^ g 
cm~^ to 1 X 10^^ g cm~^ (to 2 x 10^^ g cm~^ in some cases) in 200 logarithmically equal steps. 
This traces a mass-radius relation for a certain equation of state. We then repeat this procedure 
for a range of a to see the effect of the higher order terms in perturbative f{R) gravity. 

Recently, the authors of [43] showed that the measurement of masses and radii of three neutron 
stars are sufficient for constraining the pressure of nuclear matter at densities a few times the density 
of nuclear saturation. These data are provided by the measurements on the neutron stars EXO 
1745-248 [44], 4U 1608-52 [45] and 4U 1820-30 [46] by the methods described in the cited papers. 
We use the constraints on the M-R relation of neutron stars given in [38] , which is a union of these 
three constraints.^ The constraint of [38] is shown in all M-R plots as the region bounded by the 
thin black line. Note that the M-R constraint in the first version of [38] is larger in the published 
version. In the earlier versions of our work we employed the earlier tighter constraint and reached 
somewhat different conclusions for EPS and SLy EoS. 

Apart from the above constraint, the recent measurement [41] of the mass of the neutron star 
PSR J1614-2230 with 1.97 it 0.04 provides a stringent constraint on any M-R relation that can 
be obtained with a combination of a and EoS. This constraint is shown as the horizontal black line 
with its error shown in grey. Any viable combination of a and EoS must yield a M-R relation with 
a maximum mass exceeding this measured mass. 

The constraints on the M-R relation obtained by [38] and the 2 solar mass neutron star 
PSR J1614-2230 exclude many of the possible EoS' if one assumes GR as the ultimate classical 
theory of gravity. In the gravity model employed here, the value of a provides a new degree of 
freedom such that some of the EoS', which are excluded within the framework of GR, can now be 
reconciled with the observations for certain values of a. In the following we discuss this for all EoS' 
individually. To save space in the figures, we define the parameter ag = a/10^ cm^. We show the 
stable configurations [dM/dpc > 0) with solid lines and the unstable configurations with dashed 
lines of the same color. 

In figures, we show our results (mass versus radius, M-R, relations) for 6 representative EoS' 
for the /(i?) = R + aR^ gravity model. Results for each EoS are summarized as follows: 

• FPS (Figure 1): For FPS [47], the maximum mass that a neutron star can have, within 
GR, is about 1.8 Mq and is less than the measured mass of PSR J1614-2230. This means FPS 
can not represent the EoS of neutron stars in GR (a = 0). The maximum mass increases 
with decreasing value of a and we find that, for og < —2, the maximum mass becomes 
-^max = 2.04Mq. We thus find that FPS is consistent with the measurement of the maximum 
mass for a < —2 x 10^ cm^. Nevertheless, this does not mean a = — 10 x 10^ cm^ will be 
satisfying both constraints. In fact for such large values of |q| we can not obtain M-R relations 

^For a critic of these constraints see [39]. 
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resembling known properties of neutron stars. Furthermore for \a\ > 10^^ cm^ the vahdity of 
the perturbative approach is dubious, as we mention in the discussion section. 

• SLy (Figure 2): SLy [48] is consistent with both constraints within the framework of GR. 
For a > 2 X 10^ cm^, however, we see that Mmax is less than the measured mass of PSR 
J1614-2230. We conclude for the gravity model employed here, f{R) = R + aR^ , that SLy is 
consistent with the observations only if a < 2 x 10^ cm^. 

• AP4 (Figure 3): AP4 [49] is consistent with the constraints as long as a < 4 x 10^ cm^. 
Interestingly, we find that a new stable solution branch, for which dM/dpc > is satisfied, 
for values of a different than zero. This stable branch is obtained, for ctg = — 2 starting from 
central densities 8.6 x 10^^ g cm~^. 

• GSl (Figure 4): For GSl [50], the maximum mass in GR remains well below the measured 
mass of PSR J1614-2230. The maximum mass of neutron stars for this EoS can reach up to 
~ 2 Mq for ag = —4. Starting from ag = —2 the stability condition (dM/dpc > 0) is satisfied 
for the whole range of central densities considered. 

• MPAl (Figure 5): MPAl [51] provides a maximum mass above the observed mass of PSR 
J1614-2230 in GR, though it does not pass through the M-R constraint of [38]. For ag > 6 it 
can not satisfy the maximum mass constraint as well. 

• MSI (Figure 6): The maximum mass for MSI [52] satisfies the observed mass of PSR 
J1614-2230 only for ag < 2 though it moves away from the M-R constraint of [38] for such 
low values of a. 

For all EoS' we observe that the maximum stable mass of a neutron star, Mmax) and its radius 
at this mass, Rmin: increases for decreasing values of a, for the ranges we consider in the figures. 
There is no change in the behavior of Afmax and i?min values while a changes sign. Thus the 
structure of neutron stars in GR (a = 0) does not constitute an extremal configuration in terms of 
Mmax and Rmin- In figure 7 we show the dependence of these quantities on the value of a for the 
polytropic EoS 



used in [37] where T = 9/5 is the polytropic index and K is a constant. We fit the numerical results 
with cubic polynomials. For the maximum mass fitted with 



where Mq is the maximum mass obtained for general relativity, we find that A = —1.30796 x 10 ^ it 
5.547xlO~*^M0, B = 1.44851 x 10~'^±1.625 x 10^^ Mq and C = -3.82907x 10~3±1.234x 10"^ Mq. 
Interestingly, we find that for a = 15 x 10^ cm^ the maximum mass and minimum radius of the 
neutron star attains their minimum values starting to increase (again) beyond this value. The 
analysis is complicated by increased numerical oscillations with increasing values of |a|. Similarly, 
we fit the minimum value of NS radius which is attained at the maximum mass with a cubic function 




(3.1) 



max — 



Aog + Bal + Cag + Mq 



(3.2) 



R. 



oag + bag 



+ cag + Mq. 



(3.3) 
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We find that a = 6.18915 x 10"^ ± 1.625 x 10"*^ km, b = 0.00144416 ± 4.76 x 10^^ km, c = 
-0.0469634 ± 0.0003614 km and Rq = 11.2939 ± 0.003055 km. 

The upper and lower bounds on the value of a presented for each EoS are in the range of 
\a\ ~ 10^ cm^. Values of |a|, that are an order of magnitude smaUer than this value, produce 
results that can not be distinguished from the results obtained within GR. This corresponds to a 
curvature scale of Rq ~ a^^ ~ 10"^*^ cm~^ and a corresponding length scale of L ~ a^/^ ~ 10^ cm, 
which is an order of magnitude smaller than the radius of the neutron star. 
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Figure 1. M-R relation obtained with f{R) = i? + aR^ using the FPS. The observational constraints of 
[38] is shown with the thin black contour; the measured mass M = 1.97 ± 0.04 Afg of PSR J1614-2230 [41] is 
shown as the horizontal black line with grey errorbar. Each solid line corresponds to a stable configuration 
for a specific value of a. Dashed lines show the solutions for unstable configurations {dM/dpc < 0). The 
grey shaded region shows where the total mass would be enclosed within its Schwarzschild radius. The red 
line {a = 0) shows the result for GR. Mmax and i?min increase for decreasing values of a. Variations in the 
M-R relation comparable to employing different EoS' can be obtained for \a\ ^ 10^ cm^. Using a < 10® cm^ 
gives M-R relations that can not be distinguished from the GR result on this plot. 
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Figure 2. M-R relation for the SLy. The notation in the figure is the same as that of Figure 1 and the 
resuhs are discussed in the text. 




5 10 15 20 



R. (km) 

Figure 3. M-R relation for the AP4. The notation in the figure is the same as that of Figure 1 and the 
resuhs are discussed in the text. 
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Figure 4. M-R relation for the GSl. The notation in the figure is the same as that of Figure 1 and the 
resuhs are discussed in the text. 




Figure 5. M-R relation for the MPAl. The notation in the figure is the same as that of Figure 1 and the 
results are discussed in the text. 
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Figure 6. M-R relation for the MSI. The notation in the figure is the same as that of Figure 1 and the 
resuhs are discussed in the text. 
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Figure 7. 

Equation (3 



A/max (left 



panel) and i?min (right panel) changing with a for the polytropic EoS given in 



4 Discussion 

In this work we analyze the neutron star solutions with realistic EoS' in perturbative f{R) gravity. 
Among the modified gravity theories the f{R) theories are relatively simple to handle. However, 
even for these theories, the field equations are complicated and obtaining modified TOV equations 
in a standard fashion is difficult. This difficulty is mainly due to field equations being fourth order 
unlike in the case of general relativity, which has second order field equations. In order to resolve 
this situation, we adapt a perturbative approach [37] in which the extra terms in the gravity action 
are multiplied by a 'dimensionful' parameter a. The extra terms with some appropriate value of a 
are supposed to act perturbatively and modify the results obtained in the case of general relativity. 
We present how the perturbative f{R) modifications affect the TOV equations. 

A drawback in using neutron stars for testing alternative theories of gravity has been the 
weakly constrained M-R relation. After the tight constraints obtained in [38], it seems this is no 
longer quite true. The result of [38] and the measured mass of PSR J1614-2230 [41] excludes many 
EoS' in the framework of GR. In the f{R) = R + aR^ gravity model, the value of a provides a 
new degree of freedom and we show in this paper that some of the EoS', which are excluded within 
the framework of GR, can now be reconciled with the observations for certain values of a. This 
then brings the question of degeneracy between using different EoS' and modifying gravity. In 
the gravity model studied here, variations in M-R relation comparable to that of using different 
EoS are induced for a being in the order of 10^ cm^, which we specify for each EoS. An order 
of magnitude larger values of a, which is still 5 orders of magnitude smaller than the constraint 
obtained via Gravity Probe B [42], does not produce neutron stars with observed properties. Thus, 
we argue that |a| ^ 10^'^ cm^ is a reasonable constraint independent of the EoS. We conclude that 
the presence of uncertainties in the EoS does not cloak the effect of the free parameter a on the 
results. 

For some EoS' (AP4 and GSl) a new stable solution branch {dM/dpc > 0), which does not 
exist in general relativity, is found. This solution branch, for larger values of |a| has a larger domain 
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and joins the conventional stable branch beyond some a. In this case there is no critical maximum 
mass to the neutron stars. 

One might be curious whether the perturbative approach followed in this paper holds for the 
range of a considered. For a neutron star the typical value of the Ricci curvature is calculated to be 
roughly on the order of ~ 10"^^ cm~^. Therefore in the case of /(i?) = + ai?^ model one easily 
sees that the perturbative term is 10"^ orders of magnitude smaller than the Einstein-Hilbert term 
i?, which justifies our approach. 

Although we consider 6 representative EoS' here, the above analysis repeated with other 
realistic EoS' will not alter the order of magnitude of the constraint on a for the following reason: 
the constraints we obtained implies a length scale of i?Q ^''^ ~ 10^ cm which is only an order of 
magnitude smaller than the typical radius of a neutron star, the probe used in this test. This is 
actually the length scale below which the gravity models used here should induce modifications 
on the M-R relation of an object of size 10 km. This implies that real deviations from general 
relativity should be hidden at even much smaller values of \a\. 
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